#####Load packages
#load packages
#library(tidyverse) #there is an error with tidyverse
library(tidyr)
## Warning: package 'tidyr' was built under R version 4.1.3
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.1.3
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.1.3
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(devtools)
## Warning: package 'devtools' was built under R version 4.1.3
## Loading required package: usethis
## Warning: package 'usethis' was built under R version 4.1.3
library(stats)
library(ordinal)
## Warning: package 'ordinal' was built under R version 4.1.3
##
## Attaching package: 'ordinal'
## The following object is masked from 'package:dplyr':
##
## slice
library(MASS)
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
library(ggeffects)
## Warning: package 'ggeffects' was built under R version 4.1.3
library(effects)
## Warning: package 'effects' was built under R version 4.1.3
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.1.3
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
library(scales)
## Warning: package 'scales' was built under R version 4.1.3
library(car) #for regression diagnostics
## Warning: package 'car' was built under R version 4.1.3
##
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
##
## recode
library(broom) #for tidy output
## Warning: package 'broom' was built under R version 4.1.3
library(ggfortify) #for model diagnostics
## Warning: package 'ggfortify' was built under R version 4.1.3
#library(sjPlot) #for outputs
library(knitr) #for kable
## Warning: package 'knitr' was built under R version 4.1.3
library(effects) #for partial effects plots
library(MuMIn) #for interfacing with STAN
library(emmeans) #for estimating marginal means
## Warning: package 'emmeans' was built under R version 4.1.3
##
## Attaching package: 'emmeans'
## The following object is masked from 'package:devtools':
##
## test
library(ggeffects) #for partial effects plots
library(MASS) #for glm.nb
library(MuMIn) #for AICc
library(DHARMa) #for residual diagnostics plots
## Warning: package 'DHARMa' was built under R version 4.1.3
## This is DHARMa 0.4.6. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')
library(modelr) #for auxillary modelling functions
## Warning: package 'modelr' was built under R version 4.1.3
##
## Attaching package: 'modelr'
## The following object is masked from 'package:broom':
##
## bootstrap
## The following object is masked from 'package:ggeffects':
##
## data_grid
library(performance) #for residuals diagnostics
## Warning: package 'performance' was built under R version 4.1.3
##
## Attaching package: 'performance'
## The following objects are masked from 'package:modelr':
##
## mae, mse, rmse
library(see) #for plotting residuals
## Warning: package 'see' was built under R version 4.1.3
library(patchwork) #for grids of plots
## Warning: package 'patchwork' was built under R version 4.1.3
##
## Attaching package: 'patchwork'
## The following object is masked from 'package:MASS':
##
## area
library(brms)
## Warning: package 'brms' was built under R version 4.1.3
## Loading required package: Rcpp
## Warning: package 'Rcpp' was built under R version 4.1.3
## Loading 'brms' package (version 2.18.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').
##
## Attaching package: 'brms'
## The following object is masked from 'package:MuMIn':
##
## loo
## The following object is masked from 'package:stats':
##
## ar
library(bayesplot)
## Warning: package 'bayesplot' was built under R version 4.1.3
## This is bayesplot version 1.10.0
## - Online documentation and vignettes at mc-stan.org/bayesplot
## - bayesplot theme set to bayesplot::theme_default()
## * Does _not_ affect other ggplot2 plots
## * See ?bayesplot_theme_set for details on theme setting
library(coda)
## Warning: package 'coda' was built under R version 4.1.3
library(ggmcmc)
## Warning: package 'ggmcmc' was built under R version 4.1.3
## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2
library(broom) #for tidying outputs
library(tidybayes) #for more tidying outputs
##
## Attaching package: 'tidybayes'
## The following objects are masked from 'package:brms':
##
## dstudent_t, pstudent_t, qstudent_t, rstudent_t
library(broom.mixed)#for summarising models
## Warning: package 'broom.mixed' was built under R version 4.1.3
library(knitr)
#Arranging data set
Data 2016, treatment and control communities.
Change name of the data set
data_fiji<- read.csv ("Fiji_data_notabu_2016.4.csv",
header=TRUE,
dec=".")
Transform variables into factors (District, control/treatment)
data_fiji$village<-as.factor(data_fiji$village) #village as factor
dim(data_fiji) # n=193
## [1] 193 131
COVARIATES:
AGE
#Age (nominal)
str(data_fiji$Age)
## chr [1:193] "33" "35" "63" "69" "50" "40" "32" "34" "77" "73" "46" "60" ...
ggplot(data_fiji,aes(Age))+geom_bar()
data_fiji$Age_num<-as.numeric(as.character(data_fiji$Age))
## Warning: NAs introduced by coercion
hist(data_fiji$Age_num)
MIGRANT
#Migrant status (Categorical)
str(data_fiji$Migrant)
## int [1:193] 3 1 2 1 1 1 1 3 1 1 ...
#transform in factor
data_fiji$Migrant<-as.factor(data_fiji$Migrant)
str(data_fiji$Migrant)
## Factor w/ 4 levels "1","2","3","4": 3 1 2 1 1 1 1 3 1 1 ...
#Combine Migrant levels 2,3 and 4
data_fiji$Migrant<-as.numeric(as.character( data_fiji$Migrant))
str(data_fiji$Migrant)
## num [1:193] 3 1 2 1 1 1 1 3 1 1 ...
data_fiji$migrant2 <- as.factor(ifelse(data_fiji$Migrant > 1,2,1))
data_fiji$Gender<-as.factor(data_fiji$Gender)
#Marital status
str(data_fiji$Marital_status)
## chr [1:193] "m" "m" "m" "s" "m" "s" "s" "s" "m" "m" "m" "m" "m" "s" "w" ...
data_fiji$Marital_status<-as.factor(data_fiji$Marital_status)
#Education
str(data_fiji$Education)
## chr [1:193] "Tertiary" "Secondary " "Secondary" "Primary" "Secondary" ...
data_fiji$Education<-as.factor(data_fiji$Education)
str(data_fiji$Education)
## Factor w/ 25 levels "0","Class 2",..: 25 21 18 16 18 18 18 18 16 16 ...
data_fiji$Education_num<-as.factor(ifelse(data_fiji$Education == "Primary",1,
ifelse(data_fiji$Education=="Secondary",2,
ifelse(data_fiji$Education=="Tertiary",3,0))))
data_fiji$Education_num<-as.numeric(as.character(data_fiji$Education_num))
str(data_fiji$Education_num)
## num [1:193] 3 0 2 1 2 2 2 2 1 1 ...
#Transform in Primary, secondary and tertiary
data_fiji$Education<-recode_factor(data_fiji$Education, "Secondary-F5" = "Secondary")
levels(data_fiji$Education)
## [1] "Secondary" "0" "Class 2" "class 5"
## [5] "Class 6" "class 7" "Class 7" "class 8"
## [9] "Class 8" "Form 3" "Form 4" "Form 5"
## [13] "Form 6" "Form 7" "Intermediate" "primary"
## [17] "Primary" "secondary" "Secondary-F6" "Secondary "
## [21] "Secondary - F1" "Secondary - F4" "tertiary" "Tertiary"
data_fiji$Education<-recode_factor(data_fiji$Education, "Secondary-F6" = "Secondary")
levels(data_fiji$Education)
## [1] "Secondary" "0" "Class 2" "class 5"
## [5] "Class 6" "class 7" "Class 7" "class 8"
## [9] "Class 8" "Form 3" "Form 4" "Form 5"
## [13] "Form 6" "Form 7" "Intermediate" "primary"
## [17] "Primary" "secondary" "Secondary " "Secondary - F1"
## [21] "Secondary - F4" "tertiary" "Tertiary"
data_fiji$Education<-recode_factor(data_fiji$Education, "Secondary - F1" = "Secondary")
levels(data_fiji$Education)
## [1] "Secondary" "0" "Class 2" "class 5"
## [5] "Class 6" "class 7" "Class 7" "class 8"
## [9] "Class 8" "Form 3" "Form 4" "Form 5"
## [13] "Form 6" "Form 7" "Intermediate" "primary"
## [17] "Primary" "secondary" "Secondary " "Secondary - F4"
## [21] "tertiary" "Tertiary"
data_fiji$Education<-recode_factor(data_fiji$Education, "Secondary - F4" = "Secondary")
levels(data_fiji$Education)
## [1] "Secondary" "0" "Class 2" "class 5" "Class 6"
## [6] "class 7" "Class 7" "class 8" "Class 8" "Form 3"
## [11] "Form 4" "Form 5" "Form 6" "Form 7" "Intermediate"
## [16] "primary" "Primary" "secondary" "Secondary " "tertiary"
## [21] "Tertiary"
data_fiji$Education<-recode_factor(data_fiji$Education, "primary" = "Primary")
levels(data_fiji$Education)
## [1] "Primary" "Secondary" "0" "Class 2" "class 5"
## [6] "Class 6" "class 7" "Class 7" "class 8" "Class 8"
## [11] "Form 3" "Form 4" "Form 5" "Form 6" "Form 7"
## [16] "Intermediate" "secondary" "Secondary " "tertiary" "Tertiary"
data_fiji$Education<-recode_factor(data_fiji$Education, "secondary"= "Secondary")
levels(data_fiji$Education)
## [1] "Secondary" "Primary" "0" "Class 2" "class 5"
## [6] "Class 6" "class 7" "Class 7" "class 8" "Class 8"
## [11] "Form 3" "Form 4" "Form 5" "Form 6" "Form 7"
## [16] "Intermediate" "Secondary " "tertiary" "Tertiary"
data_fiji$Education<-recode_factor(data_fiji$Education, "tertiary"= "Tertiary")
levels(data_fiji$Education)
## [1] "Tertiary" "Secondary" "Primary" "0" "Class 2"
## [6] "class 5" "Class 6" "class 7" "Class 7" "class 8"
## [11] "Class 8" "Form 3" "Form 4" "Form 5" "Form 6"
## [16] "Form 7" "Intermediate" "Secondary "
data_fiji$Education<-recode_factor(data_fiji$Education, "Form 4"= "Tertiary")
levels(data_fiji$Education)
## [1] "Tertiary" "Secondary" "Primary" "0" "Class 2"
## [6] "class 5" "Class 6" "class 7" "Class 7" "class 8"
## [11] "Class 8" "Form 3" "Form 5" "Form 6" "Form 7"
## [16] "Intermediate" "Secondary "
data_fiji$Education<-recode_factor(data_fiji$Education, "Form 3"= "Tertiary")
levels(data_fiji$Education)
## [1] "Tertiary" "Secondary" "Primary" "0" "Class 2"
## [6] "class 5" "Class 6" "class 7" "Class 7" "class 8"
## [11] "Class 8" "Form 5" "Form 6" "Form 7" "Intermediate"
## [16] "Secondary "
data_fiji$Education<-recode_factor(data_fiji$Education, "Form 5"= "Tertiary")
levels(data_fiji$Education)
## [1] "Tertiary" "Secondary" "Primary" "0" "Class 2"
## [6] "class 5" "Class 6" "class 7" "Class 7" "class 8"
## [11] "Class 8" "Form 6" "Form 7" "Intermediate" "Secondary "
data_fiji$Education<-recode_factor(data_fiji$Education, "Form 6"= "Tertiary")
levels(data_fiji$Education)
## [1] "Tertiary" "Secondary" "Primary" "0" "Class 2"
## [6] "class 5" "Class 6" "class 7" "Class 7" "class 8"
## [11] "Class 8" "Form 7" "Intermediate" "Secondary "
data_fiji$Education<-recode_factor(data_fiji$Education, "Form 7"= "Tertiary")
levels(data_fiji$Education)
## [1] "Tertiary" "Secondary" "Primary" "0" "Class 2"
## [6] "class 5" "Class 6" "class 7" "Class 7" "class 8"
## [11] "Class 8" "Intermediate" "Secondary "
data_fiji$Education<-recode_factor(data_fiji$Education, "Class 2"= "Secondary")
levels(data_fiji$Education)
## [1] "Secondary" "Tertiary" "Primary" "0" "class 5"
## [6] "Class 6" "class 7" "Class 7" "class 8" "Class 8"
## [11] "Intermediate" "Secondary "
data_fiji$Education<-recode_factor(data_fiji$Education, "Class 5"= "Secondary")
levels(data_fiji$Education)
## [1] "Secondary" "Tertiary" "Primary" "0" "class 5"
## [6] "Class 6" "class 7" "Class 7" "class 8" "Class 8"
## [11] "Intermediate" "Secondary "
data_fiji$Education<-recode_factor(data_fiji$Education, "Class 6"= "Secondary")
levels(data_fiji$Education)
## [1] "Secondary" "Tertiary" "Primary" "0" "class 5"
## [6] "class 7" "Class 7" "class 8" "Class 8" "Intermediate"
## [11] "Secondary "
data_fiji$Education<-recode_factor(data_fiji$Education, "Class 7"= "Secondary")
levels(data_fiji$Education)
## [1] "Secondary" "Tertiary" "Primary" "0" "class 5"
## [6] "class 7" "class 8" "Class 8" "Intermediate" "Secondary "
data_fiji$Education<-recode_factor(data_fiji$Education, "class 7"= "Secondary")
levels(data_fiji$Education)
## [1] "Secondary" "Tertiary" "Primary" "0" "class 5"
## [6] "class 8" "Class 8" "Intermediate" "Secondary "
data_fiji$Education<-recode_factor(data_fiji$Education, "class 5"= "Secondary")
levels(data_fiji$Education)
## [1] "Secondary" "Tertiary" "Primary" "0" "class 8"
## [6] "Class 8" "Intermediate" "Secondary "
data_fiji$Education<-recode_factor(data_fiji$Education, "class 8"= "Secondary")
levels(data_fiji$Education)
## [1] "Secondary" "Tertiary" "Primary" "0" "Class 8"
## [6] "Intermediate" "Secondary "
data_fiji$Education<-recode_factor(data_fiji$Education, "Class 8"= "Secondary")
levels(data_fiji$Education)
## [1] "Secondary" "Tertiary" "Primary" "0" "Intermediate"
## [6] "Secondary "
data_fiji$Education<-recode_factor(data_fiji$Education, "0"= "None")
levels(data_fiji$Education)
## [1] "None" "Secondary" "Tertiary" "Primary" "Intermediate"
## [6] "Secondary "
data_fiji$Education<-recode_factor(data_fiji$Education, "Intermediate"= "Primary")
levels(data_fiji$Education)
## [1] "Primary" "None" "Secondary" "Tertiary" "Secondary "
ggplot(data_fiji,aes(Education))+geom_bar()
We will combine none and primary because there are very few people with no education
data_fiji$Education_cat3<-as.factor(ifelse(data_fiji$Education == "Tertiary",3,
ifelse(data_fiji$Education=="Secondary",2,1)))
Wealth
I have calculated wealth using all assets and chose the one with loading factors higher than 0.4 (Gurney and Darling 2017. A Global Social-Ecological Systems Monitoring Framework for Coastal Fisheries Management: A Practical Monitoring Handbook)
Loading factor is how much each variable is correlated to axis 1.
calculate wealth (Material lifestyle)
#dvd
str(data_fiji$dvd)
## chr [1:193] "y" "y" "y" "n" "n" "y" "y" "n" "y" "y" "y" "n" "n" "n" "n" ...
data_fiji$dvd<-as.factor(data_fiji$dvd)
which(is.na(data_fiji$dvd))
## integer(0)
#washing
str(data_fiji$washing)
## chr [1:193] "n" "n" "n" "n" "n" "n" "n" "n" "n" "n" "n" "n" "n" "n" "n" ...
data_fiji$washing<-as.factor(data_fiji$washing)
which(is.na(data_fiji$washing))
## integer(0)
#computer
str(data_fiji$computer)
## chr [1:193] "n" "n" "n" "n" "n" "n" "n" "n" "n" "n" "n" "n" "n" "n" "n" ...
data_fiji$computer<-as.factor(data_fiji$computer)
which(is.na(data_fiji$computer))
## integer(0)
#fridge
str(data_fiji$fridge)
## chr [1:193] "n" "n" "n" "n" "n" "n" "n" "n" "y" "n" "n" "n" "n" "n" "n" ...
data_fiji$fridge<-as.factor(data_fiji$fridge)
which(is.na(data_fiji$fridge))
## integer(0)
#tv
str(data_fiji$tv)
## chr [1:193] "y" "y" "y" "n" "n" "y" "y" "y" "y" "n" "n" "n" "n" "n" "y" ...
data_fiji$tv<-as.factor(data_fiji$tv)
which(is.na(data_fiji$tv))
## integer(0)
#satellite
str(data_fiji$satellite)
## chr [1:193] "n" "n" "n" "n" "n" "n" "n" "n" "n" "n" "n" "n" "n" "n" "n" ...
data_fiji$satellite<-as.factor(data_fiji$satellite)
data_fiji<-data_fiji%>%
mutate(satellite=na_if(satellite, "no answer"))
which(is.na(data_fiji$satellite))#187 is missing
## [1] 187
#recode rest of the variables
data_fiji$dvd<-recode_factor(data_fiji$dvd,"y" ="1")
data_fiji$dvd<-recode_factor(data_fiji$dvd,"n" ="0")
data_fiji$smartphone<-recode_factor(data_fiji$smartphone,"y" ="1")
data_fiji$smartphone<-recode_factor(data_fiji$smartphone,"n" ="0")
data_fiji$washing<-recode_factor(data_fiji$washing,"y" ="1")
data_fiji$washing<-recode_factor(data_fiji$washing,"n" ="0")
data_fiji$fridge<-recode_factor(data_fiji$fridge,"y" ="1")
data_fiji$fridge<-recode_factor(data_fiji$fridge,"n" ="0")
data_fiji$tv<-recode_factor(data_fiji$tv,"y" ="1")
data_fiji$tv<-recode_factor(data_fiji$tv,"n" ="0")
data_fiji$satellite<-recode_factor(data_fiji$satellite,"y" ="1")
data_fiji$satellite<-recode_factor(data_fiji$satellite,"n" ="0")
#transform variable into numeric vars
data_fiji$dvd<-as.numeric(as.character(data_fiji$dvd))
which(is.na(data_fiji$dvd))
## integer(0)
data_fiji$washing<-as.numeric(as.character(data_fiji$washing))
which(is.na(data_fiji$washing))
## integer(0)
data_fiji$fridge<-as.numeric(as.character(data_fiji$fridge))
which(is.na(data_fiji$fridge))
## integer(0)
data_fiji$tv<-as.numeric(as.character(data_fiji$tv))
which(is.na(data_fiji$tv))
## integer(0)
data_fiji$satellite<-as.numeric(as.character(data_fiji$satellite))
which(is.na(data_fiji$satellite)) #missing 187
## [1] 187
#PCA 18 only have one missing value.
data_fiji<-data_fiji[-c(187),]
wealth_rows_pca18.2= data_fiji[, c("dvd", "washing","fridge","tv", "satellite")]
xord18.2<-prcomp(wealth_rows_pca18.2, center= TRUE, scale=TRUE)
summary(xord18.2)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5
## Standard deviation 1.6008 1.0004 0.7855 0.67895 0.59890
## Proportion of Variance 0.5125 0.2002 0.1234 0.09219 0.07174
## Cumulative Proportion 0.5125 0.7127 0.8361 0.92826 1.00000
biplot(xord18.2)
loadings18.2<-xord18.2$rotation #rotation is the matrix of variable loadings (columns are eigenvectors) in prcomp() function
loadings18.2
## PC1 PC2 PC3 PC4 PC5
## dvd -0.4456042 0.5510784 -0.10704930 0.23661393 -0.65597547
## washing -0.4359240 -0.4797544 0.32211120 0.68440724 0.08739014
## fridge -0.4779007 -0.1855376 0.52473662 -0.66155942 -0.15549067
## tv -0.4702871 0.4877233 -0.06487272 -0.02063843 0.73233979
## satellite -0.4023035 -0.4402841 -0.77796391 -0.19368470 -0.03949986
#extract pc1
myscores.2<-xord18.2$x # "x" are the coordinates of the individuals (observations) on the principal components in prcomp()
data_fiji$wealth_msl.2<-myscores.2[,"PC1"] #add this to fiji_data and add na to the empty cells
data_fiji$wealth_msl.2
## [1] -0.7388053 -0.7388053 -0.7388053 1.2056164 1.2056164 -0.7388053
## [7] -0.7388053 0.2339693 -1.8918196 0.2328417 0.2328417 1.2056164
## [13] 1.2056164 1.2056164 0.2339693 1.2056164 -0.7388053 -0.4520388
## [19] 1.2056164 1.2056164 1.2056164 -0.7388053 1.2056164 0.2339693
## [25] 1.2056164 -0.7388053 1.2056164 0.0526021 1.2056164 0.2328417
## [31] 1.2056164 0.2339693 1.2056164 1.2056164 1.2056164 1.2056164
## [37] -1.8918196 1.2056164 -0.7388053 1.2056164 1.2056164 1.2056164
## [43] 1.2056164 1.2056164 1.2056164 1.2056164 -0.7388053 1.2056164
## [49] 1.2056164 1.2056164 -0.9049950 -0.7388053 -0.7388053 1.2056164
## [55] 1.2056164 1.2056164 -3.0307839 1.2056164 1.2056164 0.0526021
## [61] 1.2056164 1.2056164 0.2328417 1.2056164 1.2056164 1.2056164
## [67] -1.8918196 0.0666521 1.2056164 0.2339693 1.2056164 0.0526021
## [73] 1.2056164 -1.8918196 0.0666521 1.2056164 1.2056164 1.2056164
## [79] 0.0666521 0.0666521 -1.8777696 -0.7388053 1.2056164 0.0666521
## [85] 1.2056164 1.2056164 1.2056164 1.2056164 1.2056164 0.2328417
## [91] 1.2056164 0.2339693 -2.0580092 1.2056164 1.2056164 1.2056164
## [97] 0.2339693 1.2056164 1.2056164 0.2339693 -1.8918196 -2.0580092
## [103] 1.2056164 1.2056164 -3.7156644 1.2056164 0.2339693 -2.0580092
## [109] 1.2056164 -4.6884391 1.2056164 1.2056164 -0.7388053 -0.9190450
## [115] 1.2056164 -1.0863622 -3.0307839 0.0526021 -0.9190450 0.2328417
## [121] 0.0526021 -0.7388053 1.2056164 -0.7388053 -1.8777696 1.2056164
## [127] 1.2056164 1.2056164 -1.8918196 1.2056164 0.2339693 -0.7388053
## [133] 1.2056164 1.2056164 1.2056164 1.2056164 1.2056164 1.2056164
## [139] 1.2056164 1.2056164 1.2056164 -1.8918196 1.2056164 0.2328417
## [145] -0.7388053 0.2339693 -4.6884391 0.2339693 -3.7156644 -1.0863622
## [151] -1.8777696 0.0666521 -4.6884391 -3.0307839 1.2056164 -0.7388053
## [157] 1.2056164 1.2056164 -3.0307839 -2.3964605 -3.0307839 -4.6884391
## [163] -4.6884391 -2.0591368 -2.0591368 -4.6884391 -0.7388053 -3.0307839
## [169] 0.2339693 -1.8918196 1.2056164 1.2056164 1.2056164 1.2056164
## [175] -1.8918196 -1.8918196 0.0526021 1.2056164 1.2056164 1.2056164
## [181] 0.0666521 0.2339693 1.2056164 1.2056164 -4.6884391 -0.9190450
## [187] 1.2056164 1.2056164 1.2056164 1.2056164 -4.6884391 1.2056164
str(data_fiji$wealth_msl.2)
## num [1:192] -0.739 -0.739 -0.739 1.206 1.206 ...
data_fiji$wealth_msl.2_sc.1<-scale( data_fiji$wealth_msl.2, center = TRUE,scale=TRUE)
data_fiji$wealth_msl.2.resc<- rescale(data_fiji$wealth_msl.2, to = c(0,1)) #rescale
data_fiji$wealth2<-data_fiji$wealth_msl.2.resc
#gender (Gender)
str(data_fiji$Gender)
## Factor w/ 2 levels "f","m": 1 2 2 2 2 2 1 1 2 2 ...
#marital status
str(data_fiji$Marital_status)
## Factor w/ 3 levels "m","s","w": 1 1 1 2 1 2 2 2 1 1 ...
#age (Age)
str(data_fiji$Age_num)
## num [1:192] 33 35 63 69 50 40 32 34 77 73 ...
#migrant status (migrant2)
str(data_fiji$migrant2)
## Factor w/ 2 levels "1","2": 2 1 2 1 1 1 1 2 1 1 ...
#education (Education)
str(data_fiji$Education)
## Factor w/ 5 levels "Primary","None",..: 4 5 3 1 3 3 3 3 1 1 ...
str(data_fiji$Education_cat3)
## Factor w/ 3 levels "1","2","3": 3 1 2 1 2 2 2 2 1 1 ...
#wealth
str(data_fiji$wealth_msl.2)
## num [1:192] -0.739 -0.739 -0.739 1.206 1.206 ...
str(data_fiji$wealth2)
## num [1:192] 0.67 0.67 0.67 1 1 ...
Standardize continuous variables
# Scale all "numerical" variables from 0 to 1
#migrant: binomial
data_fiji$wealth2_sd <- scale(data_fiji$wealth2)*0.5
data_fiji$Age_num_sd <- scale(data_fiji$Age_num)*0.5
RESPONSE VARIABLES
I am going to create a data set for procedural equity and another one for distributional equity to avoid the elimnations of nas in procedural equity resulting from nas occuring only in disributional equity and viceversa.
DISTRIBUTIONAL EQUITY
ggplot(data_fiji, aes(Distr_Eq_Impact))+geom_bar()
Eliminate NA and don’t knows (including 6), and noanswer
data_fiji_dist <- data_fiji[data_fiji$Distr_Eq_Impact!=6,]
data_fiji_dist <- data_fiji_dist[data_fiji_dist$Distr_Eq_Impact!="dontknow",]
data_fiji_dist<- data_fiji_dist[data_fiji_dist$Distr_Eq_Impact!="noanswer",]
which(is.na(data_fiji_dist$Distr_Eq_Impact))
## integer(0)
data_fiji_dist <- data_fiji_dist[data_fiji_dist$Distr_Eq_Impact!="na",]
ggplot(data_fiji_dist, aes(Distr_Eq_Impact))+geom_bar()
ggplot(data_fiji_dist, aes(Distr_Eq_Impact,fill=Gender))+geom_bar()+theme_classic()
table<-xtabs(~Distr_Eq_Impact+migrant2+Gender, data=data_fiji_dist)
ftable(table)
## Gender f m
## Distr_Eq_Impact migrant2
## 1 1 2 2
## 2 1 1
## 2 1 6 5
## 2 7 6
## 3 1 8 2
## 2 3 1
## 4 1 19 36
## 2 17 4
## 5 1 12 19
## 2 5 4
table2<-ftable(table)
dim(data_fiji_dist)#number of observations
## [1] 160 141
which(is.na(data_fiji_dist$Age_num_sd))#checking which rows contains NAs
## [1] 43
There in 1 NA in age withing data_fiji_dist.
data_fiji_dist<-data_fiji_dist[-c(43),]#eliminate row with NA
dim(data_fiji_dist)#number of observations. N=159 for Distributional equity
## [1] 159 141
str(data_fiji_dist$Distr_Eq_Impact)
## chr [1:159] "4" "5" "5" "5" "4" "4" "4" "5" "4" "4" "3" "3" "5" "4" "4" ...
data_fiji_dist$Distr_Eq_Impact_int<-as.integer(data_fiji_dist$Distr_Eq_Impact)
str(data_fiji_dist$Distr_Eq_Impact_int)
## int [1:159] 4 5 5 5 4 4 4 5 4 4 ...
PROCEDURAL EQUITY
ggplot(data_fiji, aes(Proced_Equity))+geom_bar()
Eliminate NA and don’t knows (including 6), and noanswer
data_fiji_proced <- data_fiji[data_fiji$Proced_Equity!=6,]
data_fiji_proced <- data_fiji_proced[data_fiji_proced$Proced_Equity!="dontknow",]
data_fiji_proced<- data_fiji_proced[data_fiji_proced$Proced_Equity!="noanswer",]
which(is.na(data_fiji_proced$Proced_Equity))
## integer(0)
data_fiji_proced<- data_fiji_proced[data_fiji_proced$Proced_Equity!="na",]
ggplot(data_fiji_proced, aes(Proced_Equity))+geom_bar()
dim(data_fiji_proced)#N=168 for Procedural Equity
## [1] 168 141
str(data_fiji_proced$Proced_Equity)
## chr [1:168] "4" "5" "5" "5" "5" "4" "4" "4" "4" "3" "5" "4" "4" "4" "4" ...
data_fiji_proced$Proced_Equity_int<-as.integer(data_fiji_proced$Proced_Equity)
str(data_fiji_proced$Proced_Equity_int)
## int [1:168] 4 5 5 5 5 4 4 4 4 3 ...
Reference categories: We will chose those with more responses (migrant1, females, married, secondary education). We only have to change the reference category for Education, because the other ones are already set.
data_fiji_dist<- within(data_fiji_dist, Education_cat3 <- relevel(Education_cat3, ref = 2)) #set secondary education as the reference category as it has more counts
data_fiji_dist$Distr_Eq_int<-data_fiji_dist$Distr_Eq_Impact_int
data_fiji_dist$Distr_Eq_int_cat<-as.factor(data_fiji_dist$Distr_Eq_int)
Bayesian (brms package)
The cumulative model assumes that the effect of a predictor is equal across the rating categories (e.g. the effect of marital status is the same across the levels of distributional equity). However, this assumption may not be appropriate, and marital status may impact rating categories differently. If the effects of predictors vary across rating categories, we call the resulting model to have category specific (cs) effects.
Checking proportional odds assumption with brms
(https://discourse.mc-stan.org/t/proportional-odds-assumptions/5316)
We need to investigate if each predictor has category specific effects. In order to do this we use an adjacent category model which uses the family=acat() instead of family=cumulative() (Burkner and Vuorre xx Ordinal regression models in psychology: a tutorial https://psyarxiv.com/x8swp/).
#M.de5_add<-brm(Distr_Eq_int~Gender+migrant2+Education_cat3+wealth2_sd+Age_num_sd+Marital_status+ #(1|village),data=data_fiji_dist,family=cumulative("logit"),control = list(adapt_delta = 0.99))
#saveRDS(M.de5_add,file="M.de5_add.rda")
M.de5_add<-readRDS(file="M.de5_add.rda")
prior_summary(M.de5_add)# check default priors
## prior class coef group resp dpar nlpar lb ub
## (flat) b
## (flat) b Age_num_sd
## (flat) b Education_cat31
## (flat) b Education_cat33
## (flat) b Genderm
## (flat) b Marital_statuss
## (flat) b Marital_statusw
## (flat) b migrant22
## (flat) b wealth2_sd
## student_t(3, 0, 2.5) Intercept
## student_t(3, 0, 2.5) Intercept 1
## student_t(3, 0, 2.5) Intercept 2
## student_t(3, 0, 2.5) Intercept 3
## student_t(3, 0, 2.5) Intercept 4
## student_t(3, 0, 2.5) sd 0
## student_t(3, 0, 2.5) sd village 0
## student_t(3, 0, 2.5) sd Intercept village 0
## source
## default
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## default
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## default
## (vectorized)
## (vectorized)
Add weakly informative priors.
Weakly informative priors
priors2 <-prior(normal(0, 2.5), class = "b")+
prior(normal(0, 2.5), class = "Intercept")+
prior(student_t(3,0,2), class = "sd")
Modelling cs(Education_cat3) does not work.We will model it without specific categories for education
#```{r } M.de5_acat_all_noed<-brm(bf(Distr_Eq_int~cs(Gender)+cs(migrant2)+Education_cat3+cs(wealth2_sd)+cs(Age_num_sd)+cs(Marital_status)+ (1|village)), data=data_fiji_dist, prior=priors2, sample_prior=‘yes’, family=acat(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))
#```
#saveRDS(M.de5_acat_all_noed,file="M.de5_acat_all_noed.rda")
M.de5_acat_all_noed<-readRDS(file="M.de5_acat_all_noed.rda")
loo.M.de5_acat_all_noed<-loo(M.de5_acat_all_noed)
## Warning: Found 2 observations with a pareto_k > 0.7 in model
## 'M.de5_acat_all_noed'. It is recommended to set 'moment_match = TRUE' in order
## to perform moment matching for problematic observations.
M.de5_acat_all2<-readRDS(file="M.de5_acat_all2.rda")
loo.M.de5_acat_all2<-loo(M.de5_acat_all2)
#```{r } M.de5_acat_all2<-brm(bf(Distr_Eq_int~Gender+migrant2+Education_cat3+wealth2_sd+Age_num_sd+Marital_status+ (1|village)), data=data_fiji_dist, prior=priors2, sample_prior=‘yes’, family=acat(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))
#```
#saveRDS(M.de5_acat_all2,file="M.de5_acat_all2.rda")
M.de5_acat_all2<-readRDS(file="M.de5_acat_all2.rda")
loo.M.de5_acat_all2<-loo(M.de5_acat_all2)
Additive model
#```{r additive model2, cashe=TRUE} M.de5_add_2<-brm(bf(Distr_Eq_int~Gender+migrant2+Education_cat3+wealth2_sd+Age_num_sd+Marital_status+ (1|village)), data=data_fiji_dist, prior=priors2, sample_prior=‘yes’, family=cumulative(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))
#```
#{r} saveRDS(M.de5_add_2,file="M.de5_add_2.rda") #
M.de5_add_2<-readRDS(file="M.de5_add_2.rda")
pars <- M.de5_add_2%>% get_variables()
wch <- pars %>% stringr::str_detect("^b_((?!Intercept).)*$|^sd_.*") # Using a negative look-around
g <- purrr::map(pars[wch],
~(M.de5_add_2 %>%
hypothesis(paste(.x,"=0"), class="") %>%
plot(plot=FALSE))[[1]])
names(g) <- pars[wch]
patchwork::wrap_plots(g)
# g %>% sjPlot::plot_grid()
PP_check:
yrep_char <- posterior_predict(M.de5_add_2)
yrep_int <- sapply(data.frame(yrep_char, stringsAsFactors = TRUE), as.integer)
y_int <- data_fiji_dist$Distr_Eq_int
ppc_bars(y_int, yrep_int)
Trace plots:
M.de5_add_2%>%mcmc_trace
Chain convergence:
M.de5_add_2%>% mcmc_plot(type='rhat_hist')
## Warning: Dropped 1 NAs from 'new_rhat(rhat)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
M.de5_add_2%>% mcmc_plot(type='neff_hist')
## Warning: Dropped 1 NAs from 'new_neff_ratio(ratio)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Autocorrelation
M.de5_add_2%>% mcmc_plot(type='acf_bar')
## Warning: Removed 84 rows containing missing values (`geom_bar()`).
posterior_M.de5_add_2<- as.array(M.de5_add_2)
Residual plot:
preds<-posterior_predict(M.de5_add_2,ndraws=250,summary=FALSE)
M.de5_add_2.resids<-createDHARMa(simulatedResponse = t(preds),
observedResponse = data_fiji_dist$Distr_Eq_int,
fittedPredictedResponse = apply(preds,2,median),
# 2 = we tell R to take preds and apply it tothe rows (1) or to the columns (2)
integerResponse = FALSE)
plot(M.de5_add_2.resids) #residuals look good
summary(M.de5_add_2)
## Family: cumulative
## Links: mu = logit; disc = identity
## Formula: Distr_Eq_int ~ Gender + migrant2 + Education_cat3 + wealth2_sd + Age_num_sd + Marital_status + (1 | village)
## Data: data_fiji_dist (Number of observations: 159)
## Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 5;
## total post-warmup draws = 3200
##
## Group-Level Effects:
## ~village (Number of levels: 10)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.40 0.26 0.02 1.02 1.00 2623 3148
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1] -4.31 0.58 -5.48 -3.24 1.00 3213 3217
## Intercept[2] -2.40 0.44 -3.31 -1.55 1.00 3150 3174
## Intercept[3] -1.82 0.43 -2.69 -1.02 1.00 3209 3027
## Intercept[4] 0.57 0.40 -0.22 1.37 1.00 3105 3089
## Genderm 0.09 0.37 -0.64 0.82 1.00 2960 3180
## migrant22 -0.71 0.35 -1.39 -0.03 1.00 3277 3008
## Education_cat31 -0.85 0.42 -1.66 -0.04 1.00 3208 3173
## Education_cat33 -0.99 0.40 -1.78 -0.19 1.00 3135 3093
## wealth2_sd 0.17 0.34 -0.52 0.82 1.00 3069 2976
## Age_num_sd 0.50 0.39 -0.26 1.29 1.00 3307 3174
## Marital_statuss 0.09 0.52 -0.97 1.11 1.00 3246 3172
## Marital_statusw -1.33 0.52 -2.38 -0.32 1.00 3052 3071
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc 1.00 0.00 1.00 1.00 NA NA NA
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
loo.M.de5_add_2<-loo(M.de5_add_2)
loo.M.de5_add_2
##
## Computed from 3200 by 159 log-likelihood matrix
##
## Estimate SE
## elpd_loo -212.4 9.8
## p_loo 16.3 1.3
## looic 424.9 19.6
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
Compare cumulative and acat
loo_compare(loo.M.de5_add_2,loo.M.de5_acat_all2,loo.M.de5_acat_all_noed)
## elpd_diff se_diff
## M.de5_add_2 0.0 0.0
## M.de5_acat_all2 -1.3 1.5
## M.de5_acat_all_noed -11.6 4.9
#```{r genderxmigrant model2, cashe=TRUE} M.de5_int.mig2<-brm(bf(Distr_Eq_int~Gender*migrant2+Education_cat3+wealth2_sd+Age_num_sd+Marital_status+ (1|village)), data=data_fiji_dist, prior=priors2, sample_prior=‘yes’, family=cumulative(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))
#```
#{r} saveRDS(M.de5_int.mig2,file="M.de5_int.mig2.rda") #
M.de5_int.mig2<-readRDS(file="M.de5_int.mig2.rda")
pars <- M.de5_int.mig2 %>% get_variables()
wch <- pars %>% stringr::str_detect("^b_((?!Intercept).)*$|^sd_.*") # Using a negative look-around
g <- purrr::map(pars[wch],
~(M.de5_int.mig2 %>%
hypothesis(paste(.x,"=0"), class="") %>%
plot(plot=FALSE))[[1]])
names(g) <- pars[wch]
patchwork::wrap_plots(g)
# g %>% sjPlot::plot_grid()
PP_check:
yrep_char <- posterior_predict(M.de5_int.mig2)
yrep_int <- sapply(data.frame(yrep_char, stringsAsFactors = TRUE), as.integer)
y_int <- data_fiji_dist$Distr_Eq_int
ppc_bars(y_int, yrep_int)
Trace plots:
M.de5_int.mig2%>%mcmc_trace
Chain convergence:
M.de5_int.mig2 %>% mcmc_plot(type='rhat_hist')
## Warning: Dropped 1 NAs from 'new_rhat(rhat)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
M.de5_int.mig2 %>% mcmc_plot(type='neff_hist')
## Warning: Dropped 1 NAs from 'new_neff_ratio(ratio)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Autocorrelation
M.de5_int.mig2 %>% mcmc_plot(type='acf_bar')
## Warning: Removed 84 rows containing missing values (`geom_bar()`).
posterior_M.de5_int.mig2<- as.array(M.de5_int.mig2)
mcmc_acf(posterior_M.de5_int.mig2,
pars = c("b_Intercept[1]","b_Intercept[2]","b_Intercept[3]","b_Intercept[4]", "b_Genderm", "b_migrant22", "b_Education_cat31","b_Education_cat33", "b_wealth2_sd","b_Age_num_sd","b_Marital_statuss","b_Marital_statusw"))
Residual plot:
preds<-posterior_predict(M.de5_int.mig2,ndraws=250,summary=FALSE)
M.de5_int.mig2.resids<-createDHARMa(simulatedResponse = t(preds),
observedResponse = data_fiji_dist$Distr_Eq_int,
fittedPredictedResponse = apply(preds,2,median),
# 2 = we tell R to take preds and apply it tothe rows (1) or to the columns (2)
integerResponse = FALSE)
plot(M.de5_int.mig2.resids) #residuals look good
summary(M.de5_int.mig2)
## Family: cumulative
## Links: mu = logit; disc = identity
## Formula: Distr_Eq_int ~ Gender * migrant2 + Education_cat3 + wealth2_sd + Age_num_sd + Marital_status + (1 | village)
## Data: data_fiji_dist (Number of observations: 159)
## Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 5;
## total post-warmup draws = 3200
##
## Group-Level Effects:
## ~village (Number of levels: 10)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.43 0.26 0.03 1.04 1.00 2558 2909
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1] -4.25 0.59 -5.43 -3.16 1.00 3084 3052
## Intercept[2] -2.30 0.45 -3.22 -1.44 1.00 3215 3197
## Intercept[3] -1.70 0.43 -2.58 -0.87 1.00 3205 2963
## Intercept[4] 0.73 0.42 -0.10 1.56 1.00 3339 3116
## Genderm 0.39 0.42 -0.45 1.22 1.00 3209 3131
## migrant22 -0.27 0.43 -1.13 0.59 1.00 2985 3161
## Education_cat31 -0.90 0.42 -1.71 -0.07 1.00 3102 3175
## Education_cat33 -1.06 0.39 -1.86 -0.32 1.00 2924 3183
## wealth2_sd 0.20 0.33 -0.45 0.86 1.00 3138 3091
## Age_num_sd 0.59 0.40 -0.19 1.37 1.00 3188 3120
## Marital_statuss 0.13 0.52 -0.88 1.16 1.00 3275 3153
## Marital_statusw -1.34 0.54 -2.40 -0.29 1.00 3312 3127
## Genderm:migrant22 -1.19 0.70 -2.55 0.18 1.00 3112 3299
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc 1.00 0.00 1.00 1.00 NA NA NA
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
mcmc_intervals(posterior_M.de5_int.mig2, prob = 0.8, # 80% intervals
prob_outer = 0.95, # 95%
point_est = "mean"
,pars = c("b_Genderm", "b_migrant22", "b_Education_cat31","b_Education_cat33", "b_wealth2_sd","b_Age_num_sd","b_Marital_statuss","b_Marital_statusw","b_Genderm:migrant22"))
Conditional effects:
Ploting the interaction
conditions <- data.frame(Gender = c("f","m"))
conditional_effects(M.de5_int.mig2, "migrant2", prob=0.95,
categorical=TRUE, conditions=conditions) #re_formula=NULL takes into account random effects
inter_plot<-conditional_effects(M.de5_int.mig2, "migrant2", prob=0.95,
categorical=TRUE, conditions=conditions) #re_formula=NULL takes into account random effects
plot1<- plot(inter_plot, plot = FALSE)[[1]] + ##ggplot object
theme_classic()+scale_color_manual(values = c("1"= "darkred",
"2"= "red3",
"3"= "grey61",
"4" = "steelblue2",
"5" = "royalblue4")) +
labs(x="Migrant Status")+
theme(axis.text=element_text(size=10),
axis.title=element_text(size=15),
legend.text=element_text(size=10),
legend.title=element_text(size=10))+
scale_fill_discrete(labels = c("Very unfair", "Unfair", "Neutral","Fair", "Very fair"))
plot1
inter_plot.80<-conditional_effects(M.de5_int.mig2, "migrant2", prob=0.80,
categorical=TRUE, conditions=conditions) #re_formula=NULL takes into account random effects
plot2<- plot(inter_plot.80, plot = FALSE)[[1]] + ##ggplot object
theme_classic()+scale_color_manual(values = c("1"= "darkred",
"2"= "red3",
"3"= "grey61",
"4" = "steelblue2",
"5" = "royalblue4")) +
labs(x="Migrant Status")+
theme(axis.text=element_text(size=10),
axis.title=element_text(size=15),
legend.text=element_text(size=10),
legend.title=element_text(size=10))
plot2
Change colors
plot1
plot1<- plot(inter_plot, plot = FALSE)[[1]] + ##ggplot object
theme_classic()+scale_color_manual(values = c("1"= "#CC0033",
"2"= "#FF6699",
"3"= "#CCCCCC",
"4" = "#CCCCFF",
"5" = "#6699FF")) +
labs(x="Migrant Status")+
theme(axis.text=element_text(size=10),
axis.title=element_text(size=15),
legend.text=element_text(size=10),
legend.title=element_text(size=10))
plot1
loo.M.de5_int.mig2<-loo(M.de5_int.mig2)
loo.M.de5_int.mig2
##
## Computed from 3200 by 159 log-likelihood matrix
##
## Estimate SE
## elpd_loo -211.7 10.1
## p_loo 17.4 1.5
## looic 423.4 20.2
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
#{r M.de5_int.marital2, cashe=TRUE} M.de5_int.marital3<-brm(bf(Distr_Eq_int~Gender*Marital_status+Education_cat3+wealth2_sd+Age_num_sd+migrant2+ (1|village)), data=data_fiji_dist, prior=priors2, sample_prior='yes', family=cumulative("logit"), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99)) #
#{r} saveRDS(M.de5_int.marital3,file="M.de5_int.marital3.rda") #
M.de5_int.marital3<-readRDS(file="M.de5_int.marital3.rda")
M.de5_int.marital2<-M.de5_int.marital3
pars <- M.de5_int.marital2 %>% get_variables()
wch <- pars %>% stringr::str_detect("^b_((?!Intercept).)*$|^sd_.*") # Using a negative look-around
g <- purrr::map(pars[wch],
~(M.de5_int.marital2 %>%
hypothesis(paste(.x,"=0"), class="") %>%
plot(plot=FALSE))[[1]])
names(g) <- pars[wch]
patchwork::wrap_plots(g)
# g %>% sjPlot::plot_grid()
PP_check:
yrep_char <- posterior_predict(M.de5_int.marital2)
yrep_int <- sapply(data.frame(yrep_char, stringsAsFactors = TRUE), as.integer)
y_int <- data_fiji_dist$Distr_Eq_int
ppc_bars(y_int, yrep_int)
Trace plots:
M.de5_int.marital2%>%mcmc_trace
Chain convergence:
M.de5_int.marital2%>% mcmc_plot(type='rhat_hist')
## Warning: Dropped 1 NAs from 'new_rhat(rhat)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
M.de5_int.marital2%>% mcmc_plot(type='neff_hist')
## Warning: Dropped 1 NAs from 'new_neff_ratio(ratio)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Autocorrelation
M.de5_int.marital2%>% mcmc_plot(type='acf_bar')
## Warning: Removed 84 rows containing missing values (`geom_bar()`).
posterior_M.de5_int.marital2<- as.array(M.de5_int.marital2)
Residual plot:
preds<-posterior_predict(M.de5_int.marital2,ndraws=250,summary=FALSE)
M.de5_int.marital2.resids<-createDHARMa(simulatedResponse = t(preds),
observedResponse = data_fiji_dist$Distr_Eq_int,
fittedPredictedResponse = apply(preds,2,median),
# 2 = we tell R to take preds and apply it tothe rows (1) or to the columns (2)
integerResponse = FALSE)
plot(M.de5_int.marital2.resids) #residuals look good
summary(M.de5_int.marital2)
## Family: cumulative
## Links: mu = logit; disc = identity
## Formula: Distr_Eq_int ~ Gender * Marital_status + Education_cat3 + wealth2_sd + Age_num_sd + migrant2 + (1 | village)
## Data: data_fiji_dist (Number of observations: 159)
## Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 5;
## total post-warmup draws = 3200
##
## Group-Level Effects:
## ~village (Number of levels: 10)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.41 0.26 0.03 1.01 1.00 2371 2262
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept[1] -4.29 0.62 -5.57 -3.15 1.00 3287
## Intercept[2] -2.37 0.47 -3.32 -1.47 1.00 3093
## Intercept[3] -1.79 0.45 -2.70 -0.91 1.00 3180
## Intercept[4] 0.60 0.43 -0.21 1.46 1.00 3074
## Genderm 0.10 0.42 -0.70 0.92 1.00 3175
## Marital_statuss 0.10 0.71 -1.30 1.47 1.00 2994
## Marital_statusw -1.34 0.60 -2.53 -0.18 1.00 3119
## Education_cat31 -0.74 0.40 -1.54 0.03 1.00 3178
## Education_cat33 -0.97 0.39 -1.73 -0.20 1.00 3128
## wealth2_sd 0.16 0.34 -0.51 0.83 1.00 2859
## Age_num_sd 0.46 0.39 -0.29 1.23 1.00 3055
## migrant22 -0.70 0.36 -1.39 0.00 1.00 3015
## Genderm:Marital_statuss -0.05 0.95 -1.93 1.73 1.00 2974
## Genderm:Marital_statusw 0.08 0.92 -1.69 1.86 1.00 2526
## Tail_ESS
## Intercept[1] 2831
## Intercept[2] 3057
## Intercept[3] 3137
## Intercept[4] 2962
## Genderm 3053
## Marital_statuss 3067
## Marital_statusw 2859
## Education_cat31 3023
## Education_cat33 3115
## wealth2_sd 3027
## Age_num_sd 3039
## migrant22 2907
## Genderm:Marital_statuss 3158
## Genderm:Marital_statusw 3081
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc 1.00 0.00 1.00 1.00 NA NA NA
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
loo.M.de5_int.marital2<-loo(M.de5_int.marital2)
loo.M.de5_int.marital2
##
## Computed from 3200 by 159 log-likelihood matrix
##
## Estimate SE
## elpd_loo -215.0 9.8
## p_loo 18.4 1.6
## looic 430.1 19.6
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
#```{r M.de5_int.edu2, cashe=TRUE} M.de5_int.edu2<-brm(bf(Distr_Eq_int~Gender*Education_cat3+Marital_status+wealth2_sd+Age_num_sd+migrant2+ (1|village)), data=data_fiji_dist, prior=priors2, sample_prior=‘yes’, family=cumulative(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))
#```
#saveRDS(M.de5_int.edu2,file="M.de5_int.edu2.rda")
M.de5_int.edu2<-readRDS(file="M.de5_int.edu2.rda")
pars <- M.de5_int.edu2 %>% get_variables()
wch <- pars %>% stringr::str_detect("^b_((?!Intercept).)*$|^sd_.*") # Using a negative look-around
g <- purrr::map(pars[wch],
~(M.de5_int.edu2 %>%
hypothesis(paste(.x,"=0"), class="") %>%
plot(plot=FALSE))[[1]])
names(g) <- pars[wch]
patchwork::wrap_plots(g)
# g %>% sjPlot::plot_grid()
PP_check:
yrep_char <- posterior_predict(M.de5_int.edu2)
yrep_int <- sapply(data.frame(yrep_char, stringsAsFactors = TRUE), as.integer)
y_int <- data_fiji_dist$Distr_Eq_int
ppc_bars(y_int, yrep_int)
Trace plots:
M.de5_int.edu2%>%mcmc_trace
Chain convergence:
M.de5_int.edu2%>% mcmc_plot(type='rhat_hist')
## Warning: Dropped 1 NAs from 'new_rhat(rhat)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
M.de5_int.edu2%>% mcmc_plot(type='neff_hist')
## Warning: Dropped 1 NAs from 'new_neff_ratio(ratio)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Autocorrelation
M.de5_int.edu2%>% mcmc_plot(type='acf_bar')
## Warning: Removed 84 rows containing missing values (`geom_bar()`).
posterior_M.de5_int.edu2<- as.array(M.de5_int.edu2)
Residual plot:
preds<-posterior_predict(M.de5_int.edu2,ndraws=250,summary=FALSE)
M.de5_int.edu2.resids<-createDHARMa(simulatedResponse = t(preds),
observedResponse = data_fiji_dist$Distr_Eq_int,
fittedPredictedResponse = apply(preds,2,median),
# 2 = we tell R to take preds and apply it tothe rows (1) or to the columns (2)
integerResponse = FALSE)
plot(M.de5_int.edu2.resids) #residuals look good
summary(M.de5_int.edu2)
## Family: cumulative
## Links: mu = logit; disc = identity
## Formula: Distr_Eq_int ~ Gender * Education_cat3 + Marital_status + wealth2_sd + Age_num_sd + migrant2 + (1 | village)
## Data: data_fiji_dist (Number of observations: 159)
## Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 5;
## total post-warmup draws = 3200
##
## Group-Level Effects:
## ~village (Number of levels: 10)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.40 0.27 0.02 0.98 1.00 2691 2685
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept[1] -4.28 0.59 -5.46 -3.17 1.00 3143
## Intercept[2] -2.37 0.45 -3.26 -1.50 1.00 2702
## Intercept[3] -1.79 0.43 -2.64 -0.95 1.00 2868
## Intercept[4] 0.60 0.42 -0.20 1.43 1.00 2926
## Genderm 0.09 0.45 -0.77 0.99 1.00 3296
## Education_cat31 -0.84 0.58 -1.98 0.33 1.00 3286
## Education_cat33 -0.87 0.53 -1.89 0.13 1.00 3036
## Marital_statuss 0.09 0.54 -0.94 1.15 1.00 3022
## Marital_statusw -1.30 0.54 -2.37 -0.26 1.00 3154
## wealth2_sd 0.18 0.34 -0.51 0.87 1.00 3067
## Age_num_sd 0.48 0.40 -0.31 1.26 1.00 3496
## migrant22 -0.72 0.35 -1.40 -0.03 1.00 3202
## Genderm:Education_cat31 0.20 0.77 -1.31 1.72 1.00 3180
## Genderm:Education_cat33 -0.18 0.74 -1.58 1.24 1.00 3042
## Tail_ESS
## Intercept[1] 3303
## Intercept[2] 2992
## Intercept[3] 2988
## Intercept[4] 2848
## Genderm 3217
## Education_cat31 3221
## Education_cat33 3258
## Marital_statuss 2997
## Marital_statusw 3194
## wealth2_sd 2969
## Age_num_sd 2895
## migrant22 3174
## Genderm:Education_cat31 3112
## Genderm:Education_cat33 3010
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc 1.00 0.00 1.00 1.00 NA NA NA
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
loo.M.de5_int.edu2<-loo(M.de5_int.edu2)
loo.M.de5_int.edu2
##
## Computed from 3200 by 159 log-likelihood matrix
##
## Estimate SE
## elpd_loo -215.1 9.9
## p_loo 18.5 1.5
## looic 430.2 19.7
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
#```{r M.de5_int.age2, cashe=TRUE} M.de5_int.age2<-brm(bf(Distr_Eq_int~Gender*Age_num_sd+Education_cat3+Marital_status+wealth2_sd+migrant2+ (1|village)), data=data_fiji_dist, prior=priors2, sample_prior=‘yes’, family=cumulative(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))
#```
#{r} saveRDS(M.de5_int.age2,file="M.de5_int.age2.rda") #
M.de5_int.age2<-readRDS(file="M.de5_int.age2.rda")
pars <- M.de5_int.age2 %>% get_variables()
wch <- pars %>% stringr::str_detect("^b_((?!Intercept).)*$|^sd_.*") # Using a negative look-around
g <- purrr::map(pars[wch],
~(M.de5_int.age2 %>%
hypothesis(paste(.x,"=0"), class="") %>%
plot(plot=FALSE))[[1]])
names(g) <- pars[wch]
patchwork::wrap_plots(g)
# g %>% sjPlot::plot_grid()
PP_check:
yrep_char <- posterior_predict(M.de5_int.age2)
yrep_int <- sapply(data.frame(yrep_char, stringsAsFactors = TRUE), as.integer)
y_int <- data_fiji_dist$Distr_Eq_int
ppc_bars(y_int, yrep_int)
Trace plots:
M.de5_int.age2%>%mcmc_trace
Chain convergence:
M.de5_int.age2%>% mcmc_plot(type='rhat_hist')
## Warning: Dropped 1 NAs from 'new_rhat(rhat)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
M.de5_int.age2%>% mcmc_plot(type='neff_hist')
## Warning: Dropped 1 NAs from 'new_neff_ratio(ratio)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Autocorrelation
M.de5_int.age2%>% mcmc_plot(type='acf_bar')
## Warning: Removed 84 rows containing missing values (`geom_bar()`).
posterior_M.de5_int.age2<- as.array(M.de5_int.age2)
Residual plot:
preds<-posterior_predict(M.de5_int.age2,ndraws=250,summary=FALSE)
M.de5_int.age2.resids<-createDHARMa(simulatedResponse = t(preds),
observedResponse = data_fiji_dist$Distr_Eq_int,
fittedPredictedResponse = apply(preds,2,median),
# 2 = we tell R to take preds and apply it tothe rows (1) or to the columns (2)
integerResponse = FALSE)
plot(M.de5_int.age2.resids) #residuals look good
summary(M.de5_int.age2)
## Family: cumulative
## Links: mu = logit; disc = identity
## Formula: Distr_Eq_int ~ Gender * Age_num_sd + Education_cat3 + Marital_status + wealth2_sd + migrant2 + (1 | village)
## Data: data_fiji_dist (Number of observations: 159)
## Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 5;
## total post-warmup draws = 3200
##
## Group-Level Effects:
## ~village (Number of levels: 10)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.41 0.26 0.03 1.03 1.00 2807 3058
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1] -4.31 0.58 -5.46 -3.22 1.00 2966 3173
## Intercept[2] -2.40 0.44 -3.27 -1.54 1.00 3264 3144
## Intercept[3] -1.81 0.43 -2.67 -0.99 1.00 3191 3061
## Intercept[4] 0.58 0.41 -0.23 1.40 1.00 3319 2529
## Genderm 0.09 0.37 -0.60 0.84 1.00 3289 3127
## Age_num_sd 0.45 0.57 -0.64 1.55 1.00 2917 2972
## Education_cat31 -0.85 0.42 -1.67 -0.03 1.00 2899 3215
## Education_cat33 -0.99 0.40 -1.78 -0.21 1.00 3189 3116
## Marital_statuss 0.08 0.52 -0.96 1.07 1.00 2873 3089
## Marital_statusw -1.31 0.55 -2.39 -0.23 1.00 3052 2710
## wealth2_sd 0.17 0.35 -0.51 0.85 1.00 3329 2982
## migrant22 -0.72 0.36 -1.44 -0.04 1.00 3363 3172
## Genderm:Age_num_sd 0.08 0.66 -1.21 1.39 1.00 2693 3227
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc 1.00 0.00 1.00 1.00 NA NA NA
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
loo.M.de5_int.age2<-loo(M.de5_int.age2)
loo.M.de5_int.age2
##
## Computed from 3200 by 159 log-likelihood matrix
##
## Estimate SE
## elpd_loo -213.4 9.9
## p_loo 17.2 1.4
## looic 426.7 19.8
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
#```{r M.de5_int.wealth2, cashe=TRUE} M.de5_int.wealth2<-brm(bf(Distr_Eq_int~Gender*wealth2_sd+Age_num_sd+Education_cat3+Marital_status+migrant2+ (1|village)), data=data_fiji_dist, prior=priors2, sample_prior=‘yes’, family=cumulative(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))
#```
#saveRDS(M.de5_int.wealth2,file="M.de5_int.wealth2.rda")
M.de5_int.wealth2<-readRDS(file="M.de5_int.wealth2.rda")
pars <- M.de5_int.wealth2 %>% get_variables()
wch <- pars %>% stringr::str_detect("^b_((?!Intercept).)*$|^sd_.*") # Using a negative look-around
g <- purrr::map(pars[wch],
~(M.de5_int.wealth2 %>%
hypothesis(paste(.x,"=0"), class="") %>%
plot(plot=FALSE))[[1]])
names(g) <- pars[wch]
patchwork::wrap_plots(g)
# g %>% sjPlot::plot_grid()
PP_check:
yrep_char <- posterior_predict(M.de5_int.wealth2)
yrep_int <- sapply(data.frame(yrep_char, stringsAsFactors = TRUE), as.integer)
y_int <- data_fiji_dist$Distr_Eq_int
ppc_bars(y_int, yrep_int)
Trace plots:
M.de5_int.wealth2%>%mcmc_trace
Chain convergence:
M.de5_int.wealth2%>% mcmc_plot(type='rhat_hist')
## Warning: Dropped 1 NAs from 'new_rhat(rhat)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
M.de5_int.wealth2%>% mcmc_plot(type='neff_hist')
## Warning: Dropped 1 NAs from 'new_neff_ratio(ratio)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Autocorrelation
M.de5_int.wealth2%>% mcmc_plot(type='acf_bar')
## Warning: Removed 84 rows containing missing values (`geom_bar()`).
posterior_M.de5_int.wealth2<- as.array(M.de5_int.wealth2)
Residual plot:
preds<-posterior_predict(M.de5_int.wealth2,ndraws=250,summary=FALSE)
M.de5_int.wealth2.resids<-createDHARMa(simulatedResponse = t(preds),
observedResponse = data_fiji_dist$Distr_Eq_int,
fittedPredictedResponse = apply(preds,2,median),
# 2 = we tell R to take preds and apply it tothe rows (1) or to the columns (2)
integerResponse = FALSE)
plot(M.de5_int.wealth2.resids) #residuals look good
summary(M.de5_int.wealth2)
## Family: cumulative
## Links: mu = logit; disc = identity
## Formula: Distr_Eq_int ~ Gender * wealth2_sd + Age_num_sd + Education_cat3 + Marital_status + migrant2 + (1 | village)
## Data: data_fiji_dist (Number of observations: 159)
## Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 5;
## total post-warmup draws = 3200
##
## Group-Level Effects:
## ~village (Number of levels: 10)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.37 0.25 0.02 0.97 1.00 2687 2897
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1] -4.35 0.59 -5.61 -3.24 1.00 3063 2965
## Intercept[2] -2.42 0.44 -3.30 -1.58 1.00 2962 3015
## Intercept[3] -1.84 0.42 -2.70 -1.02 1.00 2902 2824
## Intercept[4] 0.56 0.40 -0.21 1.33 1.00 2805 3181
## Genderm 0.08 0.36 -0.62 0.77 1.00 3057 3107
## wealth2_sd 0.41 0.41 -0.38 1.21 1.00 2749 2723
## Age_num_sd 0.51 0.39 -0.25 1.27 1.00 2864 3173
## Education_cat31 -0.84 0.41 -1.66 -0.06 1.00 3220 2848
## Education_cat33 -0.97 0.39 -1.76 -0.23 1.00 2974 3117
## Marital_statuss 0.06 0.53 -0.98 1.12 1.00 3056 3303
## Marital_statusw -1.37 0.53 -2.39 -0.35 1.00 2755 2987
## migrant22 -0.71 0.35 -1.39 -0.04 1.00 3083 3134
## Genderm:wealth2_sd -0.63 0.65 -1.91 0.68 1.00 3139 2822
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc 1.00 0.00 1.00 1.00 NA NA NA
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
loo.M.de5_int.wealth2<-loo(M.de5_int.wealth2)
loo.M.de5_int.wealth2
##
## Computed from 3200 by 159 log-likelihood matrix
##
## Estimate SE
## elpd_loo -212.6 9.9
## p_loo 16.6 1.4
## looic 425.3 19.7
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
Compare models
loo_compare(loo.M.de5_int.wealth2,loo.M.de5_int.age2,loo.M.de5_int.edu2,loo.M.de5_int.marital2,loo.M.de5_int.mig2,loo.M.de5_add_2 )
## elpd_diff se_diff
## M.de5_int.mig2 0.0 0.0
## M.de5_add_2 -0.8 2.0
## M.de5_int.wealth2 -0.9 2.1
## M.de5_int.age2 -1.7 2.0
## M.de5_int.marital2 -3.3 2.1
## M.de5_int.edu2 -3.4 2.1
data_fiji_proced<- within(data_fiji_proced, Education_cat3 <- relevel(Education_cat3, ref = 2)) #set secondary education as the reference category as it has more counts
str(data_fiji_proced$Proced_Equity_int)
## int [1:168] 4 5 5 5 5 4 4 4 4 3 ...
#```{r } M.pe5_acat_all_noed<-brm(bf(Proced_Equity_int~cs(Gender)+cs(migrant2)+cs(Education_cat3)+cs(wealth2_sd)+cs(Age_num_sd)+cs(Marital_status)+ (1|village)), data=data_fiji_proced, prior=priors2, sample_prior=‘yes’, family=acat(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))
#```
#saveRDS(M.pe5_acat_all_noed,file="M.pe5_acat_all_noed.rda")
M.pe5_acat_all_noed<-readRDS(file="M.pe5_acat_all_noed.rda")
M.pe5_acat_all_noed<-loo(M.pe5_acat_all_noed)
#```{r } M.pe5_acat_all2<-brm(bf(Proced_Equity_int~Gender+migrant2+Education_cat3+wealth2_sd+Age_num_sd+Marital_status+ (1|village)), data=data_fiji_proced, prior=priors2, sample_prior=‘yes’, family=acat(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))
#```
#saveRDS(M.pe5_acat_all2,file="M.pe5_acat_all2.rda")
M.pe5_acat_all2<-readRDS(file="M.pe5_acat_all2.rda")
M.pe5_acat_all2<-readRDS(file="M.pe5_acat_all2.rda")
loo.M.pe5_acat_all2<-loo(M.pe5_acat_all2)
#```{r proced additive model2, cashe=TRUE} M.pe5_add_2<-brm(bf(Proced_Equity_int~Gender+migrant2+Education_cat3+wealth2_sd+Age_num_sd+Marital_status+ (1|village)), data=data_fiji_proced, prior=priors2, sample_prior=‘yes’, family=cumulative(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))
#```
#saveRDS(M.pe5_add_2,file="M.pe5_add_2.rda")
M.pe5_add_2<-readRDS(file="M.pe5_add_2.rda")
pars <- M.pe5_add_2%>% get_variables()
wch <- pars %>% stringr::str_detect("^b_((?!Intercept).)*$|^sd_.*") # Using a negative look-around
g <- purrr::map(pars[wch],
~(M.pe5_add_2 %>%
hypothesis(paste(.x,"=0"), class="") %>%
plot(plot=FALSE))[[1]])
names(g) <- pars[wch]
patchwork::wrap_plots(g)
# g %>% sjPlot::plot_grid()
PP_check:
yrep_char <- posterior_predict(M.pe5_add_2)
yrep_int <- sapply(data.frame(yrep_char, stringsAsFactors = TRUE), as.integer)
y_int <- data_fiji_proced$Proced_Equity_int
ppc_bars(y_int, yrep_int)
Trace plots:
M.pe5_add_2%>%mcmc_trace
Chain convergence:
M.pe5_add_2%>% mcmc_plot(type='rhat_hist')
## Warning: Dropped 1 NAs from 'new_rhat(rhat)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
M.pe5_add_2%>% mcmc_plot(type='neff_hist')
## Warning: Dropped 1 NAs from 'new_neff_ratio(ratio)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Autocorrelation
M.pe5_add_2%>% mcmc_plot(type='acf_bar')
## Warning: Removed 84 rows containing missing values (`geom_bar()`).
posterior_M.pe5_add_2<- as.array(M.pe5_add_2)
Residual plot:
preds<-posterior_predict(M.pe5_add_2,ndraws=250,summary=FALSE)
M.pe5_add_2.resids<-createDHARMa(simulatedResponse = t(preds),
observedResponse = data_fiji_proced$Proced_Equity_int,
fittedPredictedResponse = apply(preds,2,median),
# 2 = we tell R to take preds and apply it tothe rows (1) or to the columns (2)
integerResponse = FALSE)
plot(M.pe5_add_2.resids) #residuals look good
summary(M.pe5_add_2)
## Family: cumulative
## Links: mu = logit; disc = identity
## Formula: Proced_Equity_int ~ Gender + migrant2 + Education_cat3 + wealth2_sd + Age_num_sd + Marital_status + (1 | village)
## Data: data_fiji_proced (Number of observations: 168)
## Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 5;
## total post-warmup draws = 3200
##
## Group-Level Effects:
## ~village (Number of levels: 10)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.42 0.27 0.02 1.02 1.00 2559 2732
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1] -3.03 0.47 -3.99 -2.16 1.00 3324 3252
## Intercept[2] -1.86 0.41 -2.70 -1.07 1.00 3356 3172
## Intercept[3] -1.31 0.40 -2.11 -0.53 1.00 3304 3188
## Intercept[4] 0.63 0.39 -0.11 1.40 1.00 3351 3258
## Genderm 0.31 0.36 -0.40 1.00 1.00 3175 2941
## migrant22 -0.52 0.32 -1.14 0.11 1.00 3435 3256
## Education_cat31 -0.24 0.40 -1.03 0.55 1.00 3234 3161
## Education_cat33 -0.91 0.37 -1.62 -0.19 1.00 3234 3128
## wealth2_sd 0.10 0.31 -0.52 0.71 1.00 2907 2595
## Age_num_sd 0.21 0.35 -0.48 0.89 1.00 3331 3130
## Marital_statuss 0.55 0.50 -0.41 1.52 1.00 3228 3070
## Marital_statusw -0.23 0.50 -1.15 0.76 1.00 3343 3064
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc 1.00 0.00 1.00 1.00 NA NA NA
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
loo.M.pe5_add_2<-loo(M.pe5_add_2)
loo.M.pe5_add_2
##
## Computed from 3200 by 168 log-likelihood matrix
##
## Estimate SE
## elpd_loo -239.0 9.2
## p_loo 16.1 1.1
## looic 478.1 18.5
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
Compare cumulative and acat
loo_compare(loo.M.pe5_add_2,loo.M.pe5_acat_all2,M.pe5_acat_all_noed)
## elpd_diff se_diff
## M.pe5_add_2 0.0 0.0
## M.pe5_acat_all2 -0.9 1.6
## M.pe5_acat_all_noed -11.8 5.6
posterior_M.pe5_add_2<- as.array(M.pe5_add_2)
mcmc_intervals(posterior_M.pe5_add_2, prob = 0.8, # 80% intervals
prob_outer = 0.95, # 95%
point_est = "mean"
,pars = c("b_Genderm", "b_migrant22", "b_Education_cat31","b_Education_cat33", "b_wealth2_sd","b_Age_num_sd","b_Marital_statuss","b_Marital_statusw"))
#```{r proced genderxmigrant model2, cashe=TRUE} M.pe5_int.mig2<-brm(bf(Proced_Equity_int~Gender*migrant2+Education_cat3+wealth2_sd+Age_num_sd+Marital_status+ (1|village)), data=data_fiji_proced, prior=priors2, sample_prior=‘yes’, family=cumulative(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))
#```
#saveRDS(M.pe5_int.mig2,file="M.pe5_int.mig2.rda")
M.pe5_int.mig2<-readRDS(file="M.pe5_int.mig2.rda")
pars <- M.pe5_int.mig2 %>% get_variables()
wch <- pars %>% stringr::str_detect("^b_((?!Intercept).)*$|^sd_.*") # Using a negative look-around
g <- purrr::map(pars[wch],
~(M.pe5_int.mig2 %>%
hypothesis(paste(.x,"=0"), class="") %>%
plot(plot=FALSE))[[1]])
names(g) <- pars[wch]
patchwork::wrap_plots(g)
# g %>% sjPlot::plot_grid()
PP_check:
yrep_char <- posterior_predict(M.pe5_int.mig2)
yrep_int <- sapply(data.frame(yrep_char, stringsAsFactors = TRUE), as.integer)
y_int <- data_fiji_proced$Proced_Equity_int
ppc_bars(y_int, yrep_int)
Trace plots:
M.pe5_int.mig2%>%mcmc_trace
Chain convergence:
M.pe5_int.mig2 %>% mcmc_plot(type='rhat_hist')
## Warning: Dropped 1 NAs from 'new_rhat(rhat)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
M.pe5_int.mig2 %>% mcmc_plot(type='neff_hist')
## Warning: Dropped 1 NAs from 'new_neff_ratio(ratio)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Autocorrelation
M.pe5_int.mig2 %>% mcmc_plot(type='acf_bar')
## Warning: Removed 84 rows containing missing values (`geom_bar()`).
Residual plot:
preds<-posterior_predict(M.pe5_int.mig2,ndraws=250,summary=FALSE)
M.pe5_int.mig2.resids<-createDHARMa(simulatedResponse = t(preds),
observedResponse = data_fiji_proced$Proced_Equity_int,
fittedPredictedResponse = apply(preds,2,median),
# 2 = we tell R to take preds and apply it tothe rows (1) or to the columns (2)
integerResponse = FALSE)
plot(M.pe5_int.mig2.resids) #residuals look good
summary(M.pe5_int.mig2)
## Family: cumulative
## Links: mu = logit; disc = identity
## Formula: Proced_Equity_int ~ Gender * migrant2 + Education_cat3 + wealth2_sd + Age_num_sd + Marital_status + (1 | village)
## Data: data_fiji_proced (Number of observations: 168)
## Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 5;
## total post-warmup draws = 3200
##
## Group-Level Effects:
## ~village (Number of levels: 10)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.41 0.27 0.02 1.03 1.00 2280 2620
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1] -2.97 0.49 -3.96 -2.03 1.00 3280 3171
## Intercept[2] -1.81 0.43 -2.68 -0.97 1.00 3294 3165
## Intercept[3] -1.25 0.41 -2.06 -0.41 1.00 3206 3249
## Intercept[4] 0.69 0.41 -0.09 1.51 1.00 3259 3175
## Genderm 0.37 0.40 -0.41 1.17 1.00 3230 2941
## migrant22 -0.45 0.39 -1.22 0.30 1.00 3496 2510
## Education_cat31 -0.16 0.41 -0.96 0.62 1.00 3043 2908
## Education_cat33 -0.88 0.37 -1.63 -0.16 1.00 3163 3091
## wealth2_sd 0.12 0.32 -0.50 0.72 1.00 2914 3030
## Age_num_sd 0.20 0.37 -0.52 0.92 1.00 3306 2933
## Marital_statuss 0.57 0.50 -0.37 1.57 1.00 3259 3256
## Marital_statusw -0.23 0.50 -1.17 0.75 1.00 3105 2972
## Genderm:migrant22 -0.17 0.67 -1.47 1.15 1.00 3210 3285
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc 1.00 0.00 1.00 1.00 NA NA NA
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
loo.M.pe5_int.mig2<-loo(M.pe5_int.mig2)
#```{r M.pe5_int.marital2, cashe=TRUE} M.pe5_int.marital3<-brm(bf(Proced_Equity_int~Gender*Marital_status+Education_cat3+wealth2_sd+Age_num_sd+migrant2+ (1|village)), data=data_fiji_proced, prior=priors2, sample_prior=‘yes’, family=cumulative(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))
#```
#{r} saveRDS(M.pe5_int.marital3,file="M.pe5_int.marital3.rda") #
M.pe5_int.marital3<-readRDS(file="M.pe5_int.marital3.rda")
M.pe5_int.marital2<-M.pe5_int.marital3
pars <- M.pe5_int.marital2 %>% get_variables()
wch <- pars %>% stringr::str_detect("^b_((?!Intercept).)*$|^sd_.*") # Using a negative look-around
g <- purrr::map(pars[wch],
~(M.pe5_int.marital2 %>%
hypothesis(paste(.x,"=0"), class="") %>%
plot(plot=FALSE))[[1]])
names(g) <- pars[wch]
patchwork::wrap_plots(g)
# g %>% sjPlot::plot_grid()
PP_check:
yrep_char <- posterior_predict(M.pe5_int.marital2)
yrep_int <- sapply(data.frame(yrep_char, stringsAsFactors = TRUE), as.integer)
y_int <- data_fiji_proced$Proced_Equity_int
ppc_bars(y_int, yrep_int)
Trace plots:
M.pe5_int.marital2%>%mcmc_trace
Chain convergence:
M.pe5_int.marital2%>% mcmc_plot(type='rhat_hist')
## Warning: Dropped 1 NAs from 'new_rhat(rhat)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
M.pe5_int.marital2%>% mcmc_plot(type='neff_hist')
## Warning: Dropped 1 NAs from 'new_neff_ratio(ratio)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Autocorrelation
M.pe5_int.marital2%>% mcmc_plot(type='acf_bar')
## Warning: Removed 84 rows containing missing values (`geom_bar()`).
Residual plot:
preds<-posterior_predict(M.pe5_int.marital2,ndraws=250,summary=FALSE)
M.pe5_int.marital2.resids<-createDHARMa(simulatedResponse = t(preds),
observedResponse = data_fiji_proced$Proced_Equity_int,
fittedPredictedResponse = apply(preds,2,median),
# 2 = we tell R to take preds and apply it tothe rows (1) or to the columns (2)
integerResponse = FALSE)
plot(M.pe5_int.marital2.resids) #residuals look good
summary(M.pe5_int.marital2)
## Family: cumulative
## Links: mu = logit; disc = identity
## Formula: Proced_Equity_int ~ Gender * Marital_status + Education_cat3 + wealth2_sd + Age_num_sd + migrant2 + (1 | village)
## Data: data_fiji_proced (Number of observations: 168)
## Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 5;
## total post-warmup draws = 3200
##
## Group-Level Effects:
## ~village (Number of levels: 10)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.40 0.27 0.02 1.04 1.00 2905 2742
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept[1] -3.09 0.49 -4.10 -2.17 1.00 3157
## Intercept[2] -1.92 0.43 -2.78 -1.10 1.00 3093
## Intercept[3] -1.37 0.42 -2.19 -0.57 1.00 3109
## Intercept[4] 0.58 0.41 -0.20 1.38 1.00 3226
## Genderm 0.14 0.40 -0.65 0.92 1.00 3327
## Marital_statuss 0.33 0.69 -1.00 1.68 1.00 3371
## Marital_statusw -0.48 0.56 -1.55 0.62 1.00 3262
## Education_cat31 -0.11 0.41 -0.94 0.70 1.00 3240
## Education_cat33 -0.91 0.36 -1.64 -0.21 1.00 3110
## wealth2_sd 0.15 0.31 -0.47 0.79 1.00 2986
## Age_num_sd 0.22 0.36 -0.48 0.94 1.00 3421
## migrant22 -0.54 0.34 -1.21 0.11 1.00 3176
## Genderm:Marital_statuss 0.50 0.92 -1.30 2.31 1.00 3429
## Genderm:Marital_statusw 1.08 1.04 -0.90 3.14 1.00 3298
## Tail_ESS
## Intercept[1] 3194
## Intercept[2] 3024
## Intercept[3] 3162
## Intercept[4] 3090
## Genderm 3243
## Marital_statuss 3064
## Marital_statusw 3047
## Education_cat31 3169
## Education_cat33 3213
## wealth2_sd 2806
## Age_num_sd 3056
## migrant22 3296
## Genderm:Marital_statuss 2980
## Genderm:Marital_statusw 3286
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc 1.00 0.00 1.00 1.00 NA NA NA
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
loo.M.pe5_int.marital2<-loo(M.pe5_int.marital2)
loo.M.pe5_int.marital2
##
## Computed from 3200 by 168 log-likelihood matrix
##
## Estimate SE
## elpd_loo -240.5 9.2
## p_loo 17.8 1.2
## looic 480.9 18.4
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
#```{r M.pe5_int.edu2, cashe=TRUE} M.pe5_int.edu2<-brm(bf(Proced_Equity_int~Gender*Education_cat3+Marital_status+wealth2_sd+Age_num_sd+migrant2+ (1|village)), data=data_fiji_proced, prior=priors2, sample_prior=‘yes’, family=cumulative(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))
#```
#saveRDS(M.pe5_int.edu2,file="M.pe5_int.edu2.rda")
M.pe5_int.edu2<-readRDS(file="M.pe5_int.edu2.rda")
pars <- M.pe5_int.edu2 %>% get_variables()
wch <- pars %>% stringr::str_detect("^b_((?!Intercept).)*$|^sd_.*") # Using a negative look-around
g <- purrr::map(pars[wch],
~(M.pe5_int.edu2%>%
hypothesis(paste(.x,"=0"), class="") %>%
plot(plot=FALSE))[[1]])
names(g) <- pars[wch]
patchwork::wrap_plots(g)
# g %>% sjPlot::plot_grid()
PP_check:
yrep_char <- posterior_predict(M.pe5_int.edu2)
yrep_int <- sapply(data.frame(yrep_char, stringsAsFactors = TRUE), as.integer)
y_int <- data_fiji_proced$Proced_Equity_int
ppc_bars(y_int, yrep_int)
Trace plots:
M.pe5_int.edu2%>%mcmc_trace
Chain convergence:
M.pe5_int.edu2%>% mcmc_plot(type='rhat_hist')
## Warning: Dropped 1 NAs from 'new_rhat(rhat)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
M.pe5_int.edu2%>% mcmc_plot(type='neff_hist')
## Warning: Dropped 1 NAs from 'new_neff_ratio(ratio)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Autocorrelation
M.pe5_int.edu2%>% mcmc_plot(type='acf_bar')
## Warning: Removed 84 rows containing missing values (`geom_bar()`).
posterior_M.pe5_int.edu2<- as.array(M.pe5_int.edu2)
Residual plot:
preds<-posterior_predict(M.pe5_int.edu2,ndraws=250,summary=FALSE)
M.pe5_int.edu2.resids<-createDHARMa(simulatedResponse = t(preds),
observedResponse = data_fiji_proced$Proced_Equity_int,
fittedPredictedResponse = apply(preds,2,median),
# 2 = we tell R to take preds and apply it tothe rows (1) or to the columns (2)
integerResponse = FALSE)
plot(M.pe5_int.edu2.resids) #residuals look good
summary(M.pe5_int.edu2)
## Family: cumulative
## Links: mu = logit; disc = identity
## Formula: Proced_Equity_int ~ Gender * Education_cat3 + Marital_status + wealth2_sd + Age_num_sd + migrant2 + (1 | village)
## Data: data_fiji_proced (Number of observations: 168)
## Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 5;
## total post-warmup draws = 3200
##
## Group-Level Effects:
## ~village (Number of levels: 10)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.40 0.26 0.02 1.00 1.00 2776 2989
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept[1] -3.05 0.48 -4.01 -2.14 1.00 3201
## Intercept[2] -1.89 0.43 -2.77 -1.05 1.00 3144
## Intercept[3] -1.33 0.42 -2.16 -0.52 1.00 3236
## Intercept[4] 0.62 0.41 -0.16 1.43 1.00 3286
## Genderm 0.24 0.44 -0.63 1.09 1.00 3284
## Education_cat31 -0.44 0.54 -1.51 0.60 1.00 3246
## Education_cat33 -0.90 0.45 -1.79 -0.02 1.00 3196
## Marital_statuss 0.53 0.50 -0.45 1.51 1.00 3352
## Marital_statusw -0.16 0.51 -1.14 0.89 1.00 3253
## wealth2_sd 0.08 0.32 -0.56 0.70 1.00 3169
## Age_num_sd 0.21 0.36 -0.48 0.92 1.00 3069
## migrant22 -0.53 0.33 -1.18 0.11 1.00 3303
## Genderm:Education_cat31 0.45 0.77 -1.06 1.96 1.00 3321
## Genderm:Education_cat33 0.01 0.70 -1.43 1.36 1.00 3329
## Tail_ESS
## Intercept[1] 2827
## Intercept[2] 3022
## Intercept[3] 3053
## Intercept[4] 3223
## Genderm 3174
## Education_cat31 3031
## Education_cat33 3016
## Marital_statuss 3131
## Marital_statusw 3044
## wealth2_sd 3171
## Age_num_sd 3078
## migrant22 3000
## Genderm:Education_cat31 3170
## Genderm:Education_cat33 3078
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc 1.00 0.00 1.00 1.00 NA NA NA
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
loo.M.pe5_int.edu2<-loo(M.pe5_int.edu2)
loo.M.pe5_int.edu2
##
## Computed from 3200 by 168 log-likelihood matrix
##
## Estimate SE
## elpd_loo -241.1 9.3
## p_loo 18.4 1.2
## looic 482.2 18.7
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
#```{r M.pe5_int.age2, cashe=TRUE} M.pe5_int.age2<-brm(bf(Proced_Equity_int~Gender*Age_num_sd+Education_cat3+Marital_status+wealth2_sd+migrant2+ (1|village)), data=data_fiji_proced, prior=priors2, sample_prior=‘yes’, family=cumulative(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))
#``
#saveRDS(M.pe5_int.age2,file="M.pe5_int.age2.rda")
M.pe5_int.age2<-readRDS(file="M.pe5_int.age2.rda")
pars <- M.pe5_int.age2 %>% get_variables()
wch <- pars %>% stringr::str_detect("^b_((?!Intercept).)*$|^sd_.*") # Using a negative look-around
g <- purrr::map(pars[wch],
~(M.pe5_int.age2 %>%
hypothesis(paste(.x,"=0"), class="") %>%
plot(plot=FALSE))[[1]])
names(g) <- pars[wch]
patchwork::wrap_plots(g)
# g %>% sjPlot::plot_grid()
PP_check:
yrep_char <- posterior_predict(M.pe5_int.age2)
yrep_int <- sapply(data.frame(yrep_char, stringsAsFactors = TRUE), as.integer)
y_int <- data_fiji_proced$Proced_Equity_int
ppc_bars(y_int, yrep_int)
Trace plots:
M.pe5_int.age2%>%mcmc_trace
Chain convergence:
M.pe5_int.age2%>% mcmc_plot(type='rhat_hist')
## Warning: Dropped 1 NAs from 'new_rhat(rhat)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
M.pe5_int.age2%>% mcmc_plot(type='neff_hist')
## Warning: Dropped 1 NAs from 'new_neff_ratio(ratio)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Autocorrelation
M.pe5_int.age2%>% mcmc_plot(type='acf_bar')
## Warning: Removed 84 rows containing missing values (`geom_bar()`).
posterior_M.pe5_int.age2<- as.array(M.pe5_int.age2)
Residual plot:
preds<-posterior_predict(M.pe5_int.age2,ndraws=250,summary=FALSE)
M.pe5_int.age2.resids<-createDHARMa(simulatedResponse = t(preds),
observedResponse = data_fiji_proced$Proced_Equity_int,
fittedPredictedResponse = apply(preds,2,median),
# 2 = we tell R to take preds and apply it tothe rows (1) or to the columns (2)
integerResponse = FALSE)
plot(M.pe5_int.age2.resids) #residuals look good
summary(M.pe5_int.age2)
## Family: cumulative
## Links: mu = logit; disc = identity
## Formula: Proced_Equity_int ~ Gender * Age_num_sd + Education_cat3 + Marital_status + wealth2_sd + migrant2 + (1 | village)
## Data: data_fiji_proced (Number of observations: 168)
## Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 5;
## total post-warmup draws = 3200
##
## Group-Level Effects:
## ~village (Number of levels: 10)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.41 0.27 0.03 1.05 1.00 2690 2872
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1] -3.04 0.48 -3.99 -2.12 1.00 3087 2954
## Intercept[2] -1.87 0.42 -2.67 -1.05 1.00 2809 2971
## Intercept[3] -1.32 0.40 -2.09 -0.51 1.00 2981 3131
## Intercept[4] 0.63 0.40 -0.12 1.38 1.00 2962 2877
## Genderm 0.30 0.36 -0.40 1.01 1.00 3196 3054
## Age_num_sd 0.26 0.50 -0.73 1.22 1.00 2974 3053
## Education_cat31 -0.24 0.41 -1.02 0.57 1.00 3342 3290
## Education_cat33 -0.91 0.36 -1.61 -0.21 1.00 3232 3017
## Marital_statuss 0.54 0.52 -0.45 1.59 1.00 3098 2827
## Marital_statusw -0.25 0.53 -1.31 0.82 1.00 2922 2937
## wealth2_sd 0.11 0.31 -0.54 0.69 1.00 3280 3204
## migrant22 -0.52 0.32 -1.17 0.11 1.00 3290 3089
## Genderm:Age_num_sd -0.08 0.64 -1.31 1.18 1.00 2900 2926
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc 1.00 0.00 1.00 1.00 NA NA NA
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
loo.M.pe5_int.age2<-loo(M.pe5_int.age2)
loo.M.pe5_int.age2
##
## Computed from 3200 by 168 log-likelihood matrix
##
## Estimate SE
## elpd_loo -240.1 9.3
## p_loo 17.1 1.1
## looic 480.2 18.6
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
#```{r M.de5_int.wealth2, cashe=TRUE} M.pe5_int.wealth2<-brm(bf(Proced_Equity_int~Gender*wealth2_sd+Age_num_sd+Education_cat3+Marital_status+migrant2+ (1|village)), data=data_fiji_proced, prior=priors2, sample_prior=‘yes’, family=cumulative(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))
#```
#saveRDS(M.pe5_int.wealth2,file="M.pe5_int.wealth2.rda")
M.pe5_int.wealth2<-readRDS(file="M.pe5_int.wealth2.rda")
pars <- M.pe5_int.wealth2 %>% get_variables()
wch <- pars %>% stringr::str_detect("^b_((?!Intercept).)*$|^sd_.*") # Using a negative look-around
g <- purrr::map(pars[wch],
~(M.pe5_int.wealth2 %>%
hypothesis(paste(.x,"=0"), class="") %>%
plot(plot=FALSE))[[1]])
names(g) <- pars[wch]
patchwork::wrap_plots(g)
# g %>% sjPlot::plot_grid()
PP_check:
yrep_char <- posterior_predict(M.pe5_int.wealth2)
yrep_int <- sapply(data.frame(yrep_char, stringsAsFactors = TRUE), as.integer)
y_int <- data_fiji_proced$Proced_Equity_int
ppc_bars(y_int, yrep_int)
Trace plots:
M.pe5_int.wealth2%>%mcmc_trace
Chain convergence:
M.pe5_int.wealth2%>% mcmc_plot(type='rhat_hist')
## Warning: Dropped 1 NAs from 'new_rhat(rhat)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
M.pe5_int.wealth2%>% mcmc_plot(type='neff_hist')
## Warning: Dropped 1 NAs from 'new_neff_ratio(ratio)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Autocorrelation
M.pe5_int.wealth2%>% mcmc_plot(type='acf_bar')
## Warning: Removed 84 rows containing missing values (`geom_bar()`).
posterior_M.pe5_int.wealth2<- as.array(M.pe5_int.wealth2)
Residual plot:
preds<-posterior_predict(M.pe5_int.wealth2,ndraws=250,summary=FALSE)
M.pe5_int.wealth2.resids<-createDHARMa(simulatedResponse = t(preds),
observedResponse = data_fiji_proced$Proced_Equity_int,
fittedPredictedResponse = apply(preds,2,median),
# 2 = we tell R to take preds and apply it tothe rows (1) or to the columns (2)
integerResponse = FALSE)
plot(M.pe5_int.wealth2.resids) #residuals look good
summary(M.pe5_int.wealth2)
## Family: cumulative
## Links: mu = logit; disc = identity
## Formula: Proced_Equity_int ~ Gender * wealth2_sd + Age_num_sd + Education_cat3 + Marital_status + migrant2 + (1 | village)
## Data: data_fiji_proced (Number of observations: 168)
## Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 5;
## total post-warmup draws = 3200
##
## Group-Level Effects:
## ~village (Number of levels: 10)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.40 0.27 0.02 1.03 1.00 2777 2744
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1] -3.04 0.47 -3.97 -2.12 1.00 2865 2963
## Intercept[2] -1.87 0.42 -2.70 -1.05 1.00 3069 3021
## Intercept[3] -1.31 0.40 -2.10 -0.51 1.00 3137 2842
## Intercept[4] 0.63 0.40 -0.13 1.46 1.00 3284 3150
## Genderm 0.31 0.37 -0.41 1.04 1.00 3132 2800
## wealth2_sd 0.16 0.39 -0.61 0.94 1.00 2820 2527
## Age_num_sd 0.22 0.37 -0.51 0.95 1.00 3135 3091
## Education_cat31 -0.24 0.40 -1.03 0.57 1.00 3047 2839
## Education_cat33 -0.91 0.37 -1.65 -0.21 1.00 3074 2985
## Marital_statuss 0.55 0.52 -0.46 1.55 1.00 3234 2894
## Marital_statusw -0.23 0.49 -1.21 0.70 1.00 3185 3064
## migrant22 -0.52 0.33 -1.16 0.12 1.00 3134 3241
## Genderm:wealth2_sd -0.18 0.62 -1.40 1.03 1.00 3144 3004
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc 1.00 0.00 1.00 1.00 NA NA NA
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
loo.M.pe5_int.wealth2<-loo(M.pe5_int.wealth2)
loo.M.pe5_int.wealth2
##
## Computed from 3200 by 168 log-likelihood matrix
##
## Estimate SE
## elpd_loo -240.1 9.3
## p_loo 17.0 1.1
## looic 480.1 18.6
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
Compare models
loo_compare(loo.M.pe5_int.wealth2,loo.M.pe5_int.age2,loo.M.pe5_int.edu2,loo.M.pe5_int.marital2,loo.M.pe5_int.mig2,loo.M.pe5_add_2 )
## elpd_diff se_diff
## M.pe5_add_2 0.0 0.0
## M.pe5_int.wealth2 -1.0 0.3
## M.pe5_int.age2 -1.1 0.2
## M.pe5_int.mig2 -1.4 0.4
## M.pe5_int.marital2 -1.4 1.1
## M.pe5_int.edu2 -2.0 0.7